﻿using System;
using System.Collections.Generic;
using System.Runtime.InteropServices;
using System.Linq;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.Controls;
using System.IO;
using PlanningGIS.Util.Log;

namespace PlanningGIS.ArcGIS
{
    /// <summary>
    /// 概述：ESRI.ArcGIS.Geometry 帮助类<br />
    /// </summary>
    [ComVisible(true)]
    public static class GeometryHelper
    {
        /// <summary>
        /// 合并2个图形成1个图形
        /// </summary>
        /// <param name="geometry1">图形1</param>
        /// <param name="geometry2">图形2</param>
        /// <returns></returns>
        public static IGeometry Union(IGeometry geometry1, IGeometry geometry2)
        {
            ITopologicalOperator2 topoOperator1 = (ITopologicalOperator2)geometry1;
            topoOperator1.IsKnownSimple_2 = false;
            topoOperator1.Simplify();
            ITopologicalOperator2 topoOperator2 = (ITopologicalOperator2)geometry2;
            topoOperator2.IsKnownSimple_2 = false;
            topoOperator2.Simplify();
            IGeometry pGeometry = topoOperator1.Union((IGeometry)topoOperator2);
            ITopologicalOperator2 topoOperatorUnion = (ITopologicalOperator2)pGeometry;
            topoOperatorUnion.IsKnownSimple_2 = false;
            topoOperatorUnion.Simplify();
            topoOperator1 = null;
            topoOperator2 = null;
            return (IGeometry)topoOperatorUnion;
        }

        /// <summary>
        /// 合并多个图形成1个图形
        /// </summary>
        /// <param name="geometryCollection">几何集合。</param>
        /// <returns>合并后的几何。</returns>
        public static IGeometry Union(IGeometryCollection geometryCollection)
        {
            if (geometryCollection.GeometryCount == 0) return null;
            if (geometryCollection.GeometryCount == 1) return geometryCollection.get_Geometry(0);
            IGeometryBag geometryBag = geometryCollection as IGeometryBag;
            IGeometry geometry = geometryCollection.get_Geometry(0);
            geometryBag.SpatialReference = geometry.SpatialReference;
            ITopologicalOperator topo = null;
            switch (geometry.GeometryType)
            {
                case esriGeometryType.esriGeometryPolygon:
                    topo = new Polygon() as ITopologicalOperator;
                    break;
                case esriGeometryType.esriGeometryPolyline:
                    topo = new Polyline() as ITopologicalOperator;
                    break;
                case esriGeometryType.esriGeometryPoint:
                    topo = new Multipoint() as ITopologicalOperator;
                    break;
                default:
                    return null;
            }
            topo.ConstructUnion(geometryBag as IEnumGeometry);
            topo.Simplify();
            IGeometry g = topo as IGeometry;
            if (geometry.SpatialReference != null)
            {
                g.Project(geometry.SpatialReference);
            }
            return g;
        }

        /// <summary>
        /// 合并多个图形成1个图形
        /// </summary>
        /// <param name="geometrys"></param>
        /// <returns></returns>
        public static T Union<T>(List<T> geometrys)
        {
            if (geometrys == null || geometrys.Count == 0) return default(T);
            if (geometrys.Count == 1) return geometrys[0];
            IGeometryCollection geos = new GeometryBag() as IGeometryCollection;
            (geos as IGeometryBag).SpatialReference = (geometrys[0] as IGeometry).SpatialReference;
            foreach (IGeometry g in geometrys)
            {
                geos.AddGeometry(g);
            }
            return (T)Union(geos);
        }

        /// <summary>
        /// 求两个多边形的相交面
        /// </summary>
        /// <param name="polygon1">多边形1</param>
        /// <param name="polygon2">多边形2</param>
        /// <returns></returns>
        public static IPolygon Intersect(IPolygon polygon1, IPolygon polygon2)
        {
            if (polygon1 == null || polygon1.IsEmpty || polygon2 == null || polygon2.IsEmpty) return null;
            IRelationalOperator pRelationOperator = polygon1 as IRelationalOperator;
            if (pRelationOperator.Contains(polygon2))
                return polygon2;
            else if (pRelationOperator.Overlaps(polygon2))   
            {
                ITopologicalOperator3 pTopo = polygon1 as ITopologicalOperator3;
                polygon2 = pTopo.Intersect(polygon2, esriGeometryDimension.esriGeometry2Dimension) as IPolygon;
            }
            else
                return null;
            return polygon2;
        }

        /// <summary>
        /// 将指定的 WKT 格式字符串的值转换为等效的 ESRI.ArcGIS.Client.IGeometry。
        /// </summary>
        /// <param name="wkt">需要进行转换的 WKT 格式的字符串。</param>
        /// <returns>等效于 WKT 格式字符串的值的 ESRI.ArcGIS.Client.IGeometry。</returns>
        /// <remarks>目前只支持多边形的转换。</remarks>
        public static IGeometry ToGeometry(string wkt)
        {
            if (string.IsNullOrEmpty(wkt))
            {
                return null;
            }

            if (wkt.ToUpper().Contains("POLYGON"))
            {
                return ToPolygon(wkt);
            }
            return null;
        }

        /// <summary>
        /// 将指定的 WKT 格式字符串的值转换为等效的 ESRI.ArcGIS.Geometry.IPolygon。
        /// </summary>
        /// <param name="wkt">需要进行转换的 WKT 格式的字符串。</param>
        /// <returns>等效于 WKT 格式字符串的值的 ESRI.ArcGIS.Geometry.IPolygon。</returns>
        public static IPolygon ToPolygon(string wkt)
        {
            if (string.IsNullOrEmpty(wkt))
            {
                return null;
            }
            if (wkt.ToUpper().Contains("MULTIPOLYGON"))
            {
                IEnumerable<string> polygons = wkt.ToUpper().Replace("MULTIPOLYGON", "").Replace("(((", "").Replace(")))", "").Split(new string[] { ")),((" }, StringSplitOptions.RemoveEmptyEntries);
                return PolygonsToMultiPolygon(polygons);
            }
            if (wkt.ToUpper().Contains("POLYGON"))
            {
                IEnumerable<string> rings = wkt.ToUpper().Replace("POLYGON", "").Replace("((", "").Replace("))", "").Split(new string[] { "),(" }, StringSplitOptions.RemoveEmptyEntries);
                return RingsToPolygon(rings);
            }
            return null;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="polygons"></param>
        /// <returns></returns>
        private static IPolygon PolygonsToMultiPolygon(IEnumerable<string> polygons)
        {
            IGeometryCollection geometryCollection = new Polygon() as IGeometryCollection;
            foreach (string polygon in polygons)
            {
                IEnumerable<string> rings = polygon.Split(new string[] { "),(" }, StringSplitOptions.RemoveEmptyEntries);
                IPolygon pg = RingsToPolygon(rings);
                geometryCollection.AddGeometry(pg as IGeometry);
            }
            return geometryCollection as IPolygon;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="rings"></param>
        /// <returns></returns>
        private static IPolygon RingsToPolygon(IEnumerable<string> rings)
        {
            IGeometryCollection geometryCollection = new Polygon() as IGeometryCollection;
            foreach (string ring in rings)
            {
                IEnumerable<string> points = ring.Split(new char[] { ',' }, StringSplitOptions.RemoveEmptyEntries).Reverse();   // [VEPT] OGC 规范中要求外环逆时针、内环顺时针，而 ArcGIS 要求外环顺时针、内环逆时针，于是顺序要反转
                IRing r = PointsToRing(points);
                geometryCollection.AddGeometry(r as IGeometry);
            }
            return geometryCollection as IPolygon;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="points"></param>
        /// <returns></returns>
        private static IRing PointsToRing(IEnumerable<string> points)
        {
            IPointCollection pointCollection = new Ring();
            foreach (string point in points)
            {
                string[] coordinates = point.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);
                IPoint p = new Point()
                {
                    X = double.Parse(coordinates[0]),
                    Y = double.Parse(coordinates[1])
                };
                pointCollection.AddPoint(p);
            }
            return pointCollection as IRing;
        }
        
        /// <summary>
        /// 根据范围创建面空间对象
        /// </summary>
        /// <param name="xmin"></param>
        /// <param name="ymin"></param>
        /// <param name="xmax"></param>
        /// <param name="ymax"></param>
        /// <returns></returns>
        public static IPolygon ExtentToPolygon(double xmin, double ymin, double xmax, double ymax)
        {
            object missing = Type.Missing;
            IPolygon polygon = new Polygon() as IPolygon;
            IEnvelope envelop = new ESRI.ArcGIS.Geometry.Envelope() as IEnvelope;
            envelop.XMin = xmin;
            envelop.YMin = ymin;
            envelop.XMax = xmax;
            envelop.YMax = ymax;
            IPointCollection pColl = polygon as IPointCollection;
            pColl.AddPoint(envelop.LowerLeft, ref missing, ref missing);
            pColl.AddPoint(envelop.UpperLeft, ref missing, ref missing);
            pColl.AddPoint(envelop.UpperRight, ref missing, ref missing);
            pColl.AddPoint(envelop.LowerRight, ref missing, ref missing);
            pColl.AddPoint(envelop.LowerLeft, ref missing, ref missing);
            return polygon;
        }
        
        /// <summary>
        /// 根据经度获得投影文件,传入参数是经度的秒单位 仅适用于浙江测绘局 返回的是西安80 3度带的投影
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static ISpatialReference GetSpatialByX(double x)
        {
            IProjectedCoordinateSystem ps = null;
            try
            {
                ISpatialReferenceFactory2 sr = new SpatialReferenceEnvironment() as ISpatialReferenceFactory2;
                //根据图幅号进行投影
                if (x > 121.5 * 3600)
                {
                    ps = sr.CreateProjectedCoordinateSystem((int)esriSRProjCS4Type.esriSRProjCS_Xian1980_3_Degree_GK_CM_123E);
                }
                else if (x > 118.5 * 3600 && x <= 121.5 * 3600)
                {
                    ps = sr.CreateProjectedCoordinateSystem((int)esriSRProjCS4Type.esriSRProjCS_Xian1980_3_Degree_GK_CM_120E);
                    
                }
                else if (x > 115.5 * 3600 && x <= 118.5 * 3600)
                {
                    ps = sr.CreateProjectedCoordinateSystem((int)esriSRProjCS4Type.esriSRProjCS_Xian1980_3_Degree_GK_CM_117E);
                }
                return ps;
            }
            catch
            {
                return null;
            }
        }

        /// <summary>
        /// 根据经度获得投影文件,传入参数是经度的秒单位 仅适用于浙江测绘局 返回的是CGCS2000的投影
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static ISpatialReference GetCGCS2000SpatialByX(double x)
        {
            IProjectedCoordinateSystem ps = null;
            try
            {
                ISpatialReferenceFactory2 sr = new SpatialReferenceEnvironment() as ISpatialReferenceFactory2;
                //根据图幅号进行投影
                if (x > 121.5)
                {
                    ps = sr.CreateProjectedCoordinateSystem(4550);
                }
                else if (x > 118.5 && x <= 121.5)
                {
                    ps = sr.CreateProjectedCoordinateSystem(4549);

                }
                else if (x > 115.5  && x <= 118.5 )
                {
                    ps = sr.CreateProjectedCoordinateSystem(4548);
                }
                return ps;
            }
            catch
            {
                return null;
            }
        }

        /// <summary>
        /// 创建空间参考
        /// </summary>
        /// <param name="prj">空间参考信息：WKID编码；或prj文件；或prj字符串信息</param>
        /// <returns></returns>
        public static ISpatialReference CreateSpatialReference(string prj)
        {
            //构建空间参考工厂，只能通过反射构建，这是个静态的实体对象
            Type factoryType = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
            System.Object obj = Activator.CreateInstance(factoryType);
            ISpatialReferenceFactory2 fac = obj as ISpatialReferenceFactory2;
            if (prj.Trim().ToLower().EndsWith(".prj") && File.Exists(prj))
                return fac.CreateESRISpatialReferenceFromPRJFile(prj);
            else
            {
                int wkid = 0;
                if (int.TryParse(prj, out wkid))
                {
                    return fac.CreateSpatialReference(wkid);
                }
                else
                {
                    return fac.CreateESRISpatialReferenceFromPRJ(prj);
                }
            }
        }

        /// <summary>
        /// 坐标投影转换
        /// </summary>
        /// <param name="mapcontrol">地图控件</param>
        /// <param name="toprj">要转换到的投影</param>
        /// <param name="fromprj">原始的投影</param>
        /// <param name="param">转换参数，7参数</param>
        public static void TransformSR(IMapControl2 mapcontrol, string toprj, string fromprj = "", string param = "")
        {
            if (string.IsNullOrEmpty(toprj) || string.IsNullOrWhiteSpace(toprj)) return;
            ISpatialReference toRef = CreateSpatialReference(toprj);
            if (!string.IsNullOrEmpty(param) && !string.IsNullOrWhiteSpace(param)) // 要做不同椭球体转换，参数不能为空
            {
                ISpatialReference srIn = mapcontrol.SpatialReference;
                if (!string.IsNullOrEmpty(fromprj) && !string.IsNullOrWhiteSpace(fromprj))
                {
                    srIn = CreateSpatialReference(fromprj);
                }
                string[] ps = param.Split(',');
                if (ps.Length == 7)
                {
                    ICoordinateFrameTransformation pCFT = new CoordinateFrameTransformation() as ICoordinateFrameTransformation;
                    pCFT.PutParameters(double.Parse(ps[0]), double.Parse(ps[1]), double.Parse(ps[2]),
                        double.Parse(ps[3]), double.Parse(ps[4]), double.Parse(ps[5]), double.Parse(ps[6]));
                    pCFT.PutSpatialReferences(srIn, toRef);
                    pCFT.Name = "SR" + DateTime.Now.ToString("yyyyMMdd");
                    IMapGeographicTransformations mapgeotran = mapcontrol.Map as IMapGeographicTransformations;
                    IGeoTransformationOperationSet pGTSet = mapgeotran.GeographicTransformations;
                    pGTSet.Set(esriTransformDirection.esriTransformForward, pCFT);
                    pGTSet.Set(esriTransformDirection.esriTransformReverse, pCFT);
                    mapcontrol.Refresh();
                }
            }
            else
            {
                mapcontrol.SpatialReference = toRef;
                mapcontrol.Refresh();
            }
        }

        /// <summary>
        /// 得到要素的交点 专用于一个线段与一个折线的交点
        /// </summary>
        /// <param name="segment"></param>
        /// <param name="polyline"></param>
        /// <returns></returns>
        public static IPoint GetIntersection(ISegment segment, IPolyline polyline)
        {
            ISegmentCollection segmentCol = (ISegmentCollection)polyline;
            object Missing = Type.Missing;
            for (int i = 0; i <= segmentCol.SegmentCount - 1; i++)
            {
                IConstructMultipoint constructMultipoint = (IConstructMultipoint)new Multipoint();
                constructMultipoint.ConstructIntersection(segment, esriSegmentExtension.esriExtendEmbedded, segmentCol.get_Segment(i), esriSegmentExtension.esriNoExtension, out Missing, out Missing, out Missing);
                IPointCollection ptcol = (IPointCollection)constructMultipoint;
                if (ptcol.PointCount > 0)
                {
                    return ptcol.get_Point(0);
                }
            }
            return null;
        }
       
        /// <summary>
        /// 得到要素的交点  在图纸坐标中
        /// </summary>
        /// <param name="segment1"></param>
        /// <param name="segment2"></param>
        /// <param name="type1"></param>
        /// <param name="type2"></param>
        /// <returns></returns>
        public static IPoint GetIntersectionSeg(ISegment segment1, ISegment segment2, esriSegmentExtension type1, esriSegmentExtension type2)
        {
            IConstructMultipoint pConstructMultipoint = new Multipoint() as IConstructMultipoint;
            object Missing = Type.Missing;
            pConstructMultipoint.ConstructIntersection(segment1, type1, segment2, type2, out Missing, out Missing, out Missing);
            IPointCollection pTcoll = (IPointCollection)pConstructMultipoint;
            return pTcoll.get_Point(0);            
        }

        /// <summary>
        ///  二维要素转为三维
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public static IGeometry Geom2Dto3D(IGeometry geometry)
        {
            if (geometry == null) return null;
            IZAware zaware = geometry as IZAware;
            if (zaware.ZAware && zaware.ZSimple) return geometry;
            if (geometry is IPoint)
            {
                zaware.ZAware = true;
                IPoint p = geometry as IPoint;
                p.Z = 0;
                return p;
            }
            IPointCollection pc = geometry as IPointCollection;
            IZAware za = pc as IZAware;
            za.ZAware = true;
            for (int i = 0; i < pc.PointCount; i++)
            {
                IPoint p = pc.get_Point(i);
                p.Z = 0;
                pc.UpdatePoint(i, p);
            }
            return geometry;
        }

        /// <summary>
        ///  三维要素转为二维
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public static IGeometry Geom3Dto2D(IGeometry geometry)
        {
            if (geometry == null) return null;
            IZAware zaware = geometry as IZAware;
            if (!zaware.ZAware) return geometry;
            if (geometry is IPoint)
            {
                zaware.ZAware = false;
                IPoint p = geometry as IPoint;
                return p;
            }
            IPointCollection pc = geometry as IPointCollection;
            IZAware za = pc as IZAware;
            za.ZAware = false;
            for (int i = 0; i < pc.PointCount; i++)
            {
                IPoint p = pc.get_Point(i);
                pc.UpdatePoint(i, p);
            }
            return geometry;
        }

        /// <summary>
        /// 要素如果具有Z值，判断Z值是否有效，无效时对其进行剪切后赋值给新要素会出错。
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public static bool ValidateZValue(IGeometry geometry)
        {
            IPoint pnt = null;
            double dZ = 0;
            int validPntNum = 0;

            IZAware zaware = (IZAware)geometry;
            if (zaware.ZAware == true) return true;
            if (zaware.ZSimple == false)
            {
                IPointCollection pntCol = (IPointCollection)geometry;
                for (int i = 0; i < pntCol.PointCount; i++)
                {
                    pnt = pntCol.get_Point(i);
                    zaware = (IZAware)pnt;
                    if (zaware.ZSimple == true)
                    {
                        dZ = dZ + pnt.Z;
                        validPntNum += 1;
                    }
                }
                dZ = dZ / validPntNum;
                for (int j = 0; j < pntCol.PointCount; j++)
                {
                    pnt = pntCol.get_Point(j);
                    zaware = (IZAware)pnt;
                    if (zaware.ZSimple == false)
                    {
                        pnt.Z = dZ;
                        pntCol.UpdatePoint(j, pnt);
                    }
                }
            }
            return true;
        }
        
        /// <summary>
        /// 数据裁切
        /// </summary>
        /// <param name="pOriGeometry"></param>
        /// <param name="pGeomCut"></param>
        /// <param name="geomDimension"></param>
        /// <returns></returns>
        public static IGeometry CutFeatByPolygon(IGeometry pOriGeometry, IGeometry pGeomCut, esriGeometryDimension geomDimension)
        {
            if (pGeomCut == null || pOriGeometry.Dimension == esriGeometryDimension.esriGeometry0Dimension) return pOriGeometry;
            IGeometry pGeomClone = (pGeomCut as IClone).Clone() as IGeometry;
            ITopologicalOperator2 pTopoOperator = pGeomClone as ITopologicalOperator2;
            if (pTopoOperator.IsKnownSimple == false) pTopoOperator.Simplify();
            ITopologicalOperator2 pTopoOperator2 = pOriGeometry as ITopologicalOperator2;
            if (pTopoOperator2.IsKnownSimple == false) pTopoOperator2.Simplify();
            pGeomClone.Project(pOriGeometry.SpatialReference);
            IRelationalOperator pRel = pGeomClone as IRelationalOperator;
            if (pRel.Contains(pOriGeometry)) return pOriGeometry;
            IGeometry pGeom = pTopoOperator2.Intersect(pTopoOperator as IGeometry, geomDimension);
            pTopoOperator = pGeom as ITopologicalOperator2;
            pTopoOperator.IsKnownSimple_2 = false;
            pTopoOperator.Simplify();
            return pGeom;
        }

        #region 创建对象
        /// <summary>
        /// 根据X、Y创建点对象
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        public static IPoint CreatePoint(double x, double y)
        {
            return new Point()
            {
                X = x,
                Y = y
            };
        }
        
        /// <summary>
        /// 根据两点创建线对象
        /// </summary>
        /// <param name="fromx">起点X</param>
        /// <param name="fromy">起点Y</param>
        /// <param name="tox">终点X</param>
        /// <param name="toy">终点Y</param>
        /// <returns></returns>
        public static IPolyline CreatePolyline(double fromx, double fromy, double tox, double toy)
        {
            object missing = Type.Missing;
            IPointCollection ptcol = new Polyline();
            ptcol.AddPoint(new Point() { X = fromx, Y = fromy }, ref missing, ref missing);
            ptcol.AddPoint(new Point() { X = tox, Y = toy }, ref missing, ref missing);
            return ptcol as IPolyline;
        }
        #endregion


        /// <summary>
        /// 作缓冲，考虑经纬度坐标和高斯坐标
        /// </summary>
        /// <param name="geometry">要缓存的空间对象</param>
        /// <param name="buffer">缓存距离</param>
        /// <param name="spatialReference">空间参考：geometry.SpatialReference；(IFeatureClass as IGeoDataset).SpatialReference</param>
        /// <returns></returns>
        public static IGeometry GetBuffer(IGeometry geometry, double buffer, ISpatialReference spatialReference = null)
        {
            return geometry.GetBuffer(buffer, spatialReference);
        }

        /// <summary>
        /// 对线开方形缓冲
        /// </summary>
        /// <param name="polyline"></param>
        /// <param name="buffer"></param>
        /// <param name="spatialReference"></param>
        /// <returns></returns>
        public static IGeometry GetBufferSquare(IPolyline polyline, double buffer, ISpatialReference spatialReference = null)
        {
            if (spatialReference == null) spatialReference = polyline.SpatialReference;
            if (spatialReference != null)
            {
                IUnitConverter unitConverter = new UnitConverter();
                if (polyline.SpatialReference is IGeographicCoordinateSystem)
                    buffer = unitConverter.ConvertUnits(buffer, esriUnits.esriMeters, esriUnits.esriDecimalDegrees);
            }
            object Missing = Type.Missing;
            IConstructCurve constructCurve = new Polyline() as IConstructCurve;
            constructCurve.ConstructOffset(polyline, buffer / 2, ref Missing, ref Missing);
            IPolyline pPolyline1 = constructCurve as IPolyline;
            pPolyline1.SpatialReference = polyline.SpatialReference;
            ITopologicalOperator2 pTop = pPolyline1 as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            constructCurve = new Polyline() as IConstructCurve;
            constructCurve.ConstructOffset(polyline, -buffer / 2, ref Missing, ref Missing);
            IPolyline pPolyline2 = constructCurve as IPolyline;
            pPolyline2.SpatialReference = polyline.SpatialReference;
            pTop = pPolyline2 as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            IPolyline pPolyline3 = new Polyline() as IPolyline;
            pPolyline3.SpatialReference = polyline.SpatialReference;
            pPolyline3.FromPoint = pPolyline1.FromPoint;
            pPolyline3.ToPoint = pPolyline2.FromPoint;
            pTop = pPolyline3 as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            IPolyline pPolyline4 = new Polyline() as IPolyline;
            pPolyline4.SpatialReference = polyline.SpatialReference;
            pPolyline4.FromPoint = pPolyline1.ToPoint;
            pPolyline4.ToPoint = pPolyline2.ToPoint;
            pTop = pPolyline4 as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            pTop = pPolyline1 as ITopologicalOperator2;

            IGeometry pGeometry = pTop.Union(pPolyline3);
            pTop = pGeometry as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            pGeometry = pTop.Union(pPolyline2);
            pTop = pGeometry as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            pGeometry = pTop.Union(pPolyline4);
            pTop = pGeometry as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();

            IGeometryCollection pGeometries = pGeometry as IGeometryCollection;
            IPath pPath = pGeometries.get_Geometry(0) as IPath;
            IPointCollection pPnts = pPath as IPointCollection;
            IPolygon pPolygon = new Polygon() as IPolygon;
            pPolygon.SpatialReference = polyline.SpatialReference;
            (pPolygon as IPointCollection).AddPointCollection(pPnts);
            pTop = pPolygon as ITopologicalOperator2;
            pTop.IsKnownSimple_2 = false;
            pTop.Simplify();
            return pPolygon;
        }
        
        /// <summary>
        /// 线缩短
        /// </summary>
        /// <param name="polyline"></param>
        /// <param name="length">需要减去的长度</param>
        /// <param name="isRatio">所减的长度是否为与该线的长度的比</param>
        public static void Shrink(IPolyline polyline, double length, bool isRatio)
        {
            if (polyline == null || polyline.IsEmpty) return;
            IGeometryCollection geometryCln = polyline as IGeometryCollection;
            double realDistance = 0;            
            for (int k = geometryCln.GeometryCount - 1; k >= 0; k--)//因为在循环中需要移除某些不够长的线，所以要倒过来遍历
            {
                IGeometry geometry = geometryCln.get_Geometry(k);
                if (geometry.GeometryType == esriGeometryType.esriGeometryPolyline ||
                    geometry.GeometryType == esriGeometryType.esriGeometryPath)
                {
                    ICurve curve = geometry as ICurve;
                    IPointCollection pnts = geometry as IPointCollection;                    
                    if (isRatio)//先都转换成非比例尺的值
                    {
                        if (length >= 0.5)//如果要减去的部分大于等于该段的一半，则减完之后这一段就没有了，下同
                        {
                            geometryCln.RemoveGeometries(k, 1);
                            continue;
                        }
                        realDistance = curve.Length * length;
                    }
                    else
                    {
                        IUnitConverter unitConverter = new UnitConverter();
                        ISpatialReference pSpatialReference = polyline.SpatialReference;
                        if (pSpatialReference is IGeographicCoordinateSystem)
                            realDistance = unitConverter.ConvertUnits(length, esriUnits.esriMeters, esriUnits.esriDecimalDegrees);
                        if (realDistance >= curve.Length / 2)
                        {
                            geometryCln.RemoveGeometries(k, 1);
                            continue;
                        }
                    }
                    //先把新段的起点和终点确定下来，注意这里不能直接把起始点移植新的位置，否则会产生来回的折线，而是要把多余的部分先去掉
                    IPoint pointFrom = new Point();
                    curve.QueryPoint(esriSegmentExtension.esriExtendAtFrom, realDistance, false, pointFrom);
                    IPoint pointTo = new Point();
                    curve.QueryPoint(esriSegmentExtension.esriExtendAtFrom, curve.Length - realDistance, false, pointTo);
                    //去掉起点前面多余的部分并设置新起点
                    ISegmentCollection segmentCln = geometry as ISegmentCollection;
                    double len = 0;
                    int j = 0;
                    ISegment segment;
                    for (int i = 0; i < segmentCln.SegmentCount; i++)
                    {
                        segment = segmentCln.get_Segment(i);
                        len += segment.Length;
                        if (len >= realDistance) break;
                        j++;
                    }
                    if (j > 0) segmentCln.RemoveSegments(0, j, false);
                    pnts.UpdatePoint(0, pointFrom);
                    
                    //去点终点后面多余的部分并设置新终点
                    len = 0;
                    j = 0;
                    if (segmentCln.SegmentCount > 1)
                    {
                        for (int i = segmentCln.SegmentCount - 1; i >= 0; i--)
                        {
                            segment = segmentCln.get_Segment(i);
                            len += segment.Length;
                            if (len >= realDistance) break;
                            j++;
                        }
                    }
                    if (j > 0) segmentCln.RemoveSegments(segmentCln.SegmentCount - j, j, false);
                    pnts.UpdatePoint(pnts.PointCount - 1, pointTo);
                }
            }
        }

        /// <summary>
        /// 拓扑简化
        /// </summary>
        /// <param name="geometry"></param>
        public static void Simplify(IGeometry geometry)
        {
            if (geometry.GeometryType != esriGeometryType.esriGeometryPoint)
            {
                ITopologicalOperator2 top = geometry as ITopologicalOperator2;
                top.IsKnownSimple_2 = false;
                top.Simplify();
            }
        }
        
        /// <summary>
        /// 通过geometry找点
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public static IPoint GetPntFromGeom(IGeometry geometry)
        {
            if (geometry == null || geometry.IsEmpty) return null;
            IPoint pnt;
            switch (geometry.GeometryType)
            {
                case esriGeometryType.esriGeometryPoint:
                    pnt = (IPoint)geometry;
                    break;
                case esriGeometryType.esriGeometryPolyline:
                case esriGeometryType.esriGeometryPath:
                case esriGeometryType.esriGeometryLine:                    
                    pnt = new Point();
                    ICurve curve = geometry as ICurve;
                    curve.QueryPoint(esriSegmentExtension.esriNoExtension, 0.5, true, pnt);
                    break;
                case esriGeometryType.esriGeometryRing:
                case esriGeometryType.esriGeometryPolygon:
                    IArea area = (IArea)geometry;
                    try
                    {
                        pnt = area.LabelPoint;
                    }
                    catch//modified by sheng 执行pPnt = pArea.LabelPoint; 出现异常后，居然会将pGeometry的IsEmpty由False转为True，费解！
                    {                        
                        if (geometry.IsEmpty) return null;
                        pnt = area.Centroid;
                    }
                    break;
                case esriGeometryType.esriGeometryBag:
                    IGeometryCollection geomCln = geometry as IGeometryCollection;
                    pnt = GetPntFromGeom(geomCln.get_Geometry(0));
                    break;
                default:
                    pnt = null;
                    break;
            }
            return pnt;
        }

        /// <summary>
        /// 判断两个点是否为同一个点（由于两个点不可能坐标完全相等，这里用近似的方式判断两个点是否在同一位置）
        /// </summary>
        /// <param name="pnt1">点1</param>
        /// <param name="pnt2">点2</param>
        /// <returns>判断结构，相同为True；不同为False</returns>
        public static bool IsSamePnt(IPoint pnt1, IPoint pnt2)
        {
            if (Math.Abs(pnt1.X - pnt2.X) < 0.000000001)
            {
                if (Math.Abs(pnt1.Y - pnt2.Y) < 0.000000001) return true;
            }
            return false;
        }

        /// <summary>
        /// 功能：判断矩形(Min_x1, Min_y1)-(Max_x1, Max_y1)和矩形(Min_x2, Min_y2)-(Max_x2, Max_y2)是否相交。
        /// True 相交，False 不相交
        /// </summary>
        /// <param name="Min_x1"></param>
        /// <param name="Min_y1"></param>
        /// <param name="Max_x1"></param>
        /// <param name="Max_y1"></param>
        /// <param name="Min_x2"></param>
        /// <param name="Min_y2"></param>
        /// <param name="Max_x2"></param>
        /// <param name="Max_y2"></param>
        /// <returns></returns>
        public static bool IsEnvelopIntersect(double Min_x1, double Min_y1, double Max_x1, double Max_y1, double Min_x2, double Min_y2, double Max_x2, double Max_y2)
        {
            if (Min_x2 > Max_x1 || Max_x2 < Min_x1 || Min_y2 > Max_y1 || Max_y2 < Min_y1)
                return false;
            return true;
        }

        /// <summary>
        /// 判断rang1和rang2是否相交。True 相交，False 不相交
        /// </summary>
        /// <param name="rang1"></param>
        /// <param name="rang2"></param>
        /// <returns></returns>
        public static bool IsEnvelopIntersect(IEnvelope rang1, IEnvelope rang2)
        {
            return IsEnvelopIntersect(rang1.XMin, rang1.YMin, rang1.XMax, rang1.YMax, rang2.XMin, rang2.YMin, rang2.XMax, rang2.YMax);
        }

        /// <summary>
        /// 被裁切面 裁切后的面转换为不封闭线，线的起点与终点均为2面的交点
        /// </summary>
        /// <param name="polygon"></param>
        /// <param name="clipPolygon"></param>
        /// <returns></returns>
        public static IPolyline Polygon2Polyline(IPolygon polygon, IPolyline clipPolygon)
        {
            IPolyline newPolyline = new Polyline() as IPolyline;
            IProximityOperator prox = clipPolygon as IProximityOperator;
            IGeometryCollection geos = polygon as IGeometryCollection;
            for (int i = 0; i < geos.GeometryCount; i++)
            {
                IPointCollection pts = geos.get_Geometry(i) as IPointCollection;
                int idxFirst = -1, idxSecond = -1;
                for (int j = 0; j < pts.PointCount; j++)
                {
                    if (prox.ReturnDistance(pts.get_Point(j)) == 0)
                    {
                        if (idxFirst < 0) idxFirst = j;
                        else
                        {
                            if (idxSecond < 0)
                            {
                                idxSecond = j;
                                break;
                            }
                        }
                    }
                }
                IPointCollection ptsNew = new ESRI.ArcGIS.Geometry.Path();
                if (idxFirst >= 0 & idxSecond >= 0)
                {
                    for (int j = idxSecond; j < pts.PointCount - 1; j++)
                    {
                        ptsNew.AddPoint(pts.get_Point(j));
                    }
                    for (int j = 0; j <= idxFirst; j++)
                    {
                        ptsNew.AddPoint(pts.get_Point(j));
                    }
                }
                else
                {
                    ptsNew.AddPointCollection(pts);
                }
                object obj = Type.Missing;
                if (ptsNew.PointCount > 1)
                    (newPolyline as IGeometryCollection).AddGeometry(ptsNew as IGeometry, ref obj, ref obj);
            }
            return newPolyline;
        }

        /// <summary>
        /// 从面边界上裁切掉 面和线公共部分，输出结果为线（lbh）
        /// </summary>
        /// <param name="geo"></param>
        /// <param name="polyline"></param>
        /// <param name="buffer"></param>
        /// <param name="hasZ"></param>
        /// <returns></returns>
        public static IPolyline DiffFromIntersect(IGeometry geo, IPolyline polyline, double buffer, bool hasZ = false)
        {
            try
            {
                IPolyline pPolylineSrc = null;
                if (geo is IPolygon)
                {
                    ITopologicalOperator pTopo = geo as ITopologicalOperator;
                    pTopo.Simplify();
                    pPolylineSrc = pTopo.Boundary as IPolyline;
                }
                else if (geo is IPolyline)
                {
                    pPolylineSrc = geo as IPolyline;
                }
                ITopologicalOperator pTopo1 = pPolylineSrc as ITopologicalOperator;
                ITopologicalOperator pTopo2 = polyline as ITopologicalOperator;
                IPolygon pBuffer = pTopo2.Buffer(0.05) as IPolygon;
                ITopologicalOperator pTopo3 = pBuffer as ITopologicalOperator;
                IGeometry pIntersect = pTopo3.Intersect(pPolylineSrc, esriGeometryDimension.esriGeometry1Dimension);
                IGeometry pDiff = pTopo1.Difference(pIntersect);
                if (hasZ)
                    PlanningGIS.ArcGIS.GeometryHelper.Geom2Dto3D(pDiff);
                return pDiff as IPolyline;
            }
            catch (System.Exception ex)
            {
                LogHelper.WriteErrorText(ex);
                return null;
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="geo"></param>
        /// <param name="polyline"></param>
        /// <param name="hasZ"></param>
        /// <returns></returns>
        public static IPolyline DiffFromIntersect(IGeometry geo, IPolyline polyline, bool hasZ = false)
        {
            try
            {
                IPolyline pPolyline1 = null;
                if (geo is IPolygon)
                {
                    ITopologicalOperator pTopo = geo as ITopologicalOperator;
                    pTopo.Simplify();
                    pPolyline1 = pTopo.Boundary as IPolyline;
                }
                else if (geo is IPolyline)
                {
                    pPolyline1 = geo as IPolyline;
                }
                IPolyline pPolyline2 = polyline;
                ITopologicalOperator pTopo2 = pPolyline1 as ITopologicalOperator;
                IGeometry pIntersect = pTopo2.Intersect(polyline, esriGeometryDimension.esriGeometry1Dimension);
                IGeometry pDiff = pTopo2.Difference(pIntersect);
                if (hasZ)
                    PlanningGIS.ArcGIS.GeometryHelper.Geom2Dto3D(pDiff);
                return pDiff as IPolyline;
            }
            catch (System.Exception ex)
            {
                LogHelper.WriteErrorText(ex);
                return null;
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="lines"></param>
        /// <param name="spoint"></param>
        /// <param name="epoint"></param>
        /// <returns></returns>
        public static IPolyline InterceptLine(List<IPolyline> lines, IPoint spoint, IPoint epoint)
        {
            if (lines == null || lines.Count == 0) return null;
            IPolyline line = Union(lines);
            if (line == null || line.IsEmpty) return null;
            bool happened = false;
            int partIndex;
            int segmentIndex;
            line.SplitAtPoint(spoint, true, false, out happened, out partIndex, out segmentIndex);
            if (!happened) return null;
            bool happened2 = false;
            int partIndex2;
            int segmentIndex2;
            line.SplitAtPoint(epoint, true, false, out happened2, out partIndex2, out segmentIndex2);
            if (!happened2) return null;
            int si = segmentIndex;
            int ei = segmentIndex2;
            if (segmentIndex > segmentIndex2)
            {
                ei = segmentIndex + 1;
                si = segmentIndex2;
            }
            IPointCollection pts = line as IPointCollection;
            IPolyline newline = new Polyline() as IPolyline;
            IPointCollection pts2 = newline as IPointCollection;
            for (int i = si; i <= ei; i++)
            {
                pts2.AddPoint(pts.get_Point(i));
            }
            return newline;
        }

    }
}
